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Using the replica method, we develop an analytical approach to compute the characteristic func¬ 
tion for the probability Pn{K,X) that a large N x N adjacency matrix of sparse random graphs 
has K eigenvalues below a threshold A. The method allows to determine, in principle, all moments 
of Pn{K, a), from which the typical sample to sample fluctuations can be fully characterized. For 
random graph models with localized eigenvectors, we show that the index variance scales linearly 
with A 3> 1 for |A| > 0, with a model-dependent prefactor that can be exactly calculated. Explicit 
results are discussed for Erdos-Renyi and regular random graphs, both exhibiting a prefactor with 
a non-monotonic behavior as a function of A. These results contrast with rotationally invariant 
random matrices, where the index variance scales only as In A, with an universal prefactor that is 
independent of A. Numerical diagonalization results confirm the exactness of our approach and, in 
addition, strongly support the Gaussian nature of the index fluctuations. 
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I. INTRODUCTION 

Since the pioneering work of Wigner in the statistics 
of nuclear energy levels [I[, random matrix theory has 
established itself as a research field on its own, with 
many important applications in physics and beyond Q. 
Valuable information on the behavior of different systems 
may be extracted from the eigenvalue statistics of related 
random matrix models. In this respect, meaningful sta¬ 
tistical observables are the eigenvalue distribution, the 
distribution of extreme eigenvalues and the nearest-level 
spacing distribution, to name just a few 0. 

Another prominent observable is the index /Cat (A) of a 
N X N random matrix, defined here as the total number 
of eigenvalues below a threshold A. The random variable 
/Cjv(A) is of fundamental importance in the character¬ 
ization of disordered systems described by a potential 
energy surface ... ,xn) in the A-dimensional con¬ 

figurational space [3|. The eigenvalues of the symmet¬ 
ric Hessian matrix M, formed by the second derivatives 
Mij = d'^'H/dxidxj, encode all information regarding the 
stability properties. The number of positive (negative) 
eigenvalues counts the number of stable (unstable) direc¬ 
tions around a certain configuration, while the magnitude 
of an eigenvalue quantifies the surface curvature along 
the corresponding direction. In particular, the minima 
(maxima) of the potential energy are stationary points 
in which all Hessian eigenvalues are positive (negative). 
The index is a valuable tool to probe the energy land- 
sc^e of systems as diverse as liquids Si, spin-glasses 
i-i, synchronization models i and biomolecules Q. 

The simplest model for the Hessian of a disordered sys¬ 
tem consists in neglecting its dependency with respect to 
the configurations and assuming that the elements My 
are independently drawn from a Gaussian distribution. 
In this case, the Hessian belongs to the GOE ensemble 
of random matrices i and the index statistics has been 
studied originally in reference Q, using a fermionic ver¬ 
sion of the replica method. The authors have obtained 


the large-A behavior of the index distribution Vn{K, A) 

P^(A,A)^exp|-2^[/^-iVm(A)f |, (1) 

where m{X)= dX'p{X') follows from the Wigner semi¬ 
circle law [2 for the eigenvalue distribution p{X). Equa¬ 
tion m implies that, for A ^ 1, the index variance scales 
logarithmically with A and the typical fluctuations on a 
scale of width 0(-\/ln A) around the average index have 
a Gaussian form. 

Recently, a significant amount of work has been de¬ 
voted to study the index distribution of rotationally in¬ 
variant ensembles, including Gaussian Mini , Wishart 
M and Cauchy random matrices HI- These models 
share the property that the joint probability distribution 
of eigenvalues is analytically known, which allows to em¬ 
ploy the Coulomb gas technique, pioneered by Dyson 0, 
to compute not only the typical index distribution, but 
also its large deviation regime, which characterizes atypi¬ 
cal large fluctuations |10l - [l^ . For all these ensembles, eq. 
0 is recovered in the regime of small fluctuations, with a 
variance that grows as In A for large A. The prefactor 
cr^ is given by cr^ = I/tt^ for both Gaussian [3,[l3,[lll and 
Wishart (l^ random matrices, independently of A, while 
cr^ = 2/7r^ for Cauchy random matrices [l^. This loga¬ 
rithmic behavior of the variance apparently reflects the 
repulsion between neighboring levels [l5j | , which imposes 
a constraint on the total number of eigenvalues that fit 
in a finite region of the spectrum. 

Despite the success of the Coulomb gas approach, the 
analytical form of the joint probability distribution of 
eigenvalues is not known for various interesting random 
matrix models. Perhaps the most representative exam¬ 
ple in this sense is the adjacency matrix of sparse ran¬ 
dom graphs MM, in which the average total num¬ 
ber of nonzero entries scales only linearly with A. Al¬ 
though the eigenvalue distribution of random graphs has 
been computed using different techniques [l^ , the statis- 
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tical properties of the index have not been addressed so 
far. Several random graph models typically contain local¬ 
ized eigenvectors at finite sectors of the spectrum [311, 
usually corresponding to extreme eigenvalues, where the 
nearest-level spacing distribution follows a Poisson law 
(Hill- In these regions, neighboring eigenvalues are free 
to be arbitrarily close to each other, which should heav¬ 
ily influence the index fluctuations. Models in which the 
state variables are placed on the nodes of random graphs 
have found an enormous number of applications, includ¬ 
ing spin-glasses, satisfiability problems, error-correcting 
codes and complex networks (see [^, and references 
therein), and alternative tools to study their index fluc¬ 
tuations would be more than welcome. 

In this paper we derive an analytical expression for 
the characteristic function of the index distribution de¬ 
scribing the adjacency matrix of a broad class of random 
graphs, defined in terms of an arbitrary degree distribu¬ 
tion. In principle, such analytical result allows to cal¬ 
culate the leading contribution in the large-limit of 
all moments of Vn^KjX), yet we concentrate here on 
the first and second moments. Specifically, we show 
that the index variance of random graphs scales gen¬ 
erally as a'^{X)N, with a prefactor o-^(A) that depends 
on the threshold A and on the particular structure of 
the random graph model at hand. For random regular 
graphs with uniform edges, in which all eigenvectors are 
delocalized [3 - l3| . we show that cr^(A) = 0 for any A. 
On the other hand, for random graph models with lo¬ 
calized eigenvectors [3 - [3 .13. the prefactor ct^(A) 
exhibits a maximum for a certain A, while it vanishes 
for |A| —>■ 0. These results indicate that the linear scal¬ 
ing of the variance is a consequence of the uncorrelated 
nature of the eigenvalues in the localized regions of the 
spectrum. Since cr^(O) = 0 for random graphs with an 
arbitrary degree distribution, the linear scaling breaks 
down for A = 0 and the logarithmic scaling reemerges as 
the large-leading contribution for the index variance, 
which is supported by numerical diagonalization results. 
The model-dependent character of o-^(A) contrasts with 
the highly universal prefactor found in rotationally in¬ 
variant ensembles, though the typical index fluctuations 
of random graphs remain Gaussian distributed, as sup¬ 
ported by numerical diagonalization results. 

In the next section, we lay the ground for the replica 
computation of the characteristic function. The random 
graph model is introduced in section uni the replica ap¬ 
proach is developed in section HVl and the final analytical 
result for the characteristic function is presented in sec¬ 
tion El We discuss explicit results for the average and 
the variance of the index in section EH and, in the final 
section, some final remarks are presented. 


II. THE GENERAL SETTING 


terms of a calculation reminiscent from the statistical 
mechanics of disordered systems. Let us consider a N x 
N real symmetric matrix A with eigenvalues Ai,..., Xn- 
The density of eigenvalues between A' and A' -I- dX' reads 

N 

Pn{X') = Y, • (2) 

Q —1 

The index is defined here as the total number of eigen¬ 
values smaller than a threshold A 

pX N 

ICn{X)= dX'pN{X') = Y^i^->^c.), (3) 


where 0(...) is the Heavside step function. The object 
/Cat (A) is also regarded as the integrated density of states 
or the cumulative distribution function. At this point we 
introduce the generating function 


ZNiz) 



dcj) exp 




( 4 ) 


with (p = {pi ,..., pi^) and z = X — ie, where e > 0 is 
a regularizer that ensures the convergence of the above 
Gaussian integral and I denotes the identity matrix. The 
vector components pi,, pN are real-valued. By using 
an identity that relates the Heavside function with the 
complex logarithm, eq. m can be written in terms of 
Zm{z) as follows 

K.n{X) = —: lim ^-a.Ziq{z*)— \d.Zm{z)] . (5) 

m £->0+ 

Equation ([S]) holds for a single matrix A with an arbi¬ 
trary dimension N. 

An ensemble of random matrices is defined by a large 
set of instances of A drawn independently from a distri¬ 
bution p{A). In this paper, we are interested in comput¬ 
ing the averaged index distribution 


Vn{K,X) = {6[K-ICn{X)]), (6) 

where (...) denotes the ensemble average with p{A). Us¬ 
ing an integral representation of the Dirac delta and sub¬ 
stituting eq. in eq. we obtain 

Vn{K, X) = J A), (7) 

where the characteristic function 

C/Ar(/r,A)= lim ([Zn{z)]~‘^ [Zn{z*)]‘^') (8) 

contains the whole information about the statistical 
properties of the index. The moments of the index dis¬ 
tribution are determined from 


In this section we show how to recast the problem of 
computing the index distribution of a random matrix in 


{K-) = {-ty 


,d'^gN{p,x) 


dp" 


n e N. 


( 9 ) 


^—0 






3 


The aim here is to compute the leading contribution to 
C/Ar(^, A) for iV —>• oo. According to eq. dS]), A) is 

calculated from the ensemble average of a function that 
contains real powers of the generating function, which is 
an unfeasible computation. In order to proceed further, 
we invoke the main strategy of the replica method and 
rewrite eq. ([5]) as follows 

C/Ar(/x, A)= lim lim {[Zn{z)\^~' [ZN{z*y\'*) . (10) 

£^0+ n±->±* 


The idea is to treat initially n_ and n_|_ as integers, which 
allows to compute the ensemble average. Once this aver¬ 
age is calculated and the limit iV —> oo is taken, we make 
an analytical continuation of n± to the real values ±^. 


III. RANDOM GRAPHS WITH AN 
ARBITRARY DEGREE DISTRIBUTION 


We study the index distribution oi N x N symmetric 
adjacency matrices with the following entries 


Ajj Cij Jij , 


( 11 ) 


where Cij = Cji and = Jji. The variables Cij € {0,1} 
encode the topology of the underlying random graph: we 
set Cij = 1 if there is an edge between nodes i and j, and 
zero otherwise. The real variable Jij denotes the weight 
or the strength of the undirected coupling between the 
adjacent nodes i and j. 

Both types of random variables are drawn indepen¬ 
dently from probability distributions. At this stage, there 
is no need to specify the distribution P{J) of the entries 
Jij and the model definitions are kept as general as pos¬ 
sible. However, we do need to specify the distribution of 
{cij}, which is given by [3l| 




X 



■ N 

j=i 



where the product Y[i<j runs over all distinct pairs of 
nodes and Cat is the normalization factor. 

In this model, the topology of the corresponding graph 
is solely determined by the degree fci({ci<j })=Ef=iG, 
of each node i, defined as the total number of edges at¬ 
tached to i. According to eq. d, any two nodes are 
connected with probability c/N, in which c is the aver¬ 
age degree, while the term involving the Kronecker delta 
ensures that the number of edges attached to a certain 
node i is constrained to an integer ki. For A —>■ oo, aver¬ 
aged quantities with respect to p({ci<j}) should depend 
only upon the degree distribution 


Pk 


1 ^ 

2=1 


(13) 


Equation (I12p comprises a large class of random graph 
models with distinct degree distributions, provided they 
fulfill c = 'Yl’^=oPkk- Although the ensemble average in 
the replica approach is performed with the distribution of 
eq. dm) and the final expression for Gn{p, A) is presented 
in its full generality, we discuss in section explicit 
results for regular and Erdos-Renyi (ER) random graphs, 
where the degree distributions are given, respectively, by 

Pk = 4,c [II and Pk = HI]. 


IV. THE REPLICA APPROACH 

According to eq. (iia, the characteristic function is 
obtained by calculating the moments of the generating 
function. Substituting eq. (SD in eq. m , we can rewrite 

Nn_ Nn^ 

S„(^,A)= mn ■ (^) ■ V„{z), 

( 14 ) 


in which we have defined the function 

N 


= J J-({</),,i/)J), 


(15) 


with 


A44^)=exp + 

'»/’*}) = ( exp i^CijJij (4.0J 
\ \ i<j 



The objects 4 = (</)),..., </)” ) and ,..., 

are the replicated vectors at node i. The ensemble 
average (...) includes the average over the distribu¬ 
tion of {cij}, defined in eq. (na, and the average 
over the weights {Jij}, whose distribution P{J) is ar¬ 
bitrary. In this section we evaluate the leading term of 
N~^ InVn^iz) for A —>• oo by means of the saddle-point 
method. 

Using an integral representation for the Kronecker 
delta in eq. (El), the average over the topological dis¬ 
order is explicitly calculated and the function F reads 


\ ii=l 

where 

A(0, Ip] <p', Ip') = (exp [iJ {cp.(p' - ^p.■^p')])J , (17) 
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and ■)j stands for the average over J. We have re¬ 
tained only the leading contribution of 0{N) in the expo¬ 
nent of eq. (HU). To proceed further, the order-parameter 

1 ^ 

' i=l 


is introduced in eq. (fTOl) by means of a functional delta, 
yielding the expression 




Nc 

e 2 


C 


N 


VpVp 


X exp J dijidip p{4>,xj})p{4>-,'^) 
f cN f 

xexp(— / dcfjdtp p{(f>,xp)r{cl),tp) 


N 


dxi 


N 


n?f 


Akn 


exp|^-i^e 


(19) 


with 

r{ 4 >, '^)= J dcf>'dil}'A{4>, Ip] cp\ xp')p{cp\ xp'). (20) 


The conjugated order parameter p{(p,xp) has been 
rescaled according to p{(p, xp) —>■ Np{(p, xp) and the func¬ 
tional measure in the above integral may be written as 
VpVp = ■^dp{(p,xp)dp{cp,xp), where the product 

runs over all possible values of cp and xp. By substituting 
the large-leading contribution to Ctv in eq. (HU) 


Cn = exp 


N 



c- y^Pfclnfc! 
fc=o / 


0 ( 1 ) 


( 21 ) 


where the order-parameters p{(p,xp) and p{(p,xp) fulfill 
the saddle-point equations 

p{(p,xp) = icr{(p,xp), (25) 

p{cp xp)=y^ fcpfc H^{(p,xp)[r{(p,xp)f~^ 

fc=o J dcp'dxp'H:,{(p',xp')[r{(p',xp')]’" 

(26) 

Equations (|25|l and (|26)) are obtained by extremizing the 
action S[p, p] with respect to p and p, respectively. In¬ 
serting eqs. (1251) and (l26ll back into eq. (l23ll and noting 
from eq. (l26ll that 


dcpdxp p{cp,xp)r{cp,xp) = 1, 


we derive the compact expression 


S[pJ] = 

k^O 


d(pdxpH^{(p,xp) [r{(p,xp)] 



(27) 

The last step consists in performing the limit n± —>■ 
in the above equation. In order to make progress in this 
task, we need to make an assumption regarding the struc¬ 
ture of p{cp, xp) in the replica space. 


V. THE CHARACTERISTIC FUNCTION OF 
THE INDEX DISTRIBUTION 


We follow previous works [2^, [s^ and, with a mod¬ 
est amount of foresight, we assume that p{cp,xp) has the 
following Gaussian form 


and then inserting the resulting expression into eq. (Ha, , , 

we arrive at the integral form p{4>j '4’) — 


'Dn±{z) = / VpVpexp{NS[p,p]), 


( 22 ) 


U{n±) 

< (^) 

\27ru J 


du dv Wn±. {u, v) 


\27TU j 


I ^ ^2 ^ / 2 


(28) 


where the action reads 

5'[p,/3] = I-clnc-bi J dtpdxp p{cp,xp)p{cp,xp) 

+ ^J dcpdxp p{(p,xp)r{(p,xp) 

CO 

+X! 

k=0 

The integral in eq. (1^^ can be suitably evaluated 
through the saddle-point method. In the limit A —>■ oo, 
the function {z) is given by 


dcpdxp H^{cp,xp) [-ip{cp,xp)] 


(23) 


(24) 


where Wn^ (u, v) is the normalized joint distribution of 
the complex variances u and v, with Imit > 0 and 
Imr; > 0. The latter conditions ensure the convergence 
of the integrals in eq. (l28l) . Since p{(p,xp) is not normal¬ 
ized for arbitrary n± (see eq. (US])), the factor U{n±) 
has been consistently included in eq. (PSI) . The above 
replica symmetric (RS) form of p{(p,xp) remains invari¬ 
ant under rotations of the vectors cp and xp as well as 
under permutations of the vector components. A rigor¬ 
ous approach [s^ for the eigenvalue distribution of sparse 
random graphs has confirmed the exactness of the results 
obtained via the RS assumption. 

By inserting eq. (1^51) in eq. (1251) and then taking the 
limit n± —> ±^, one derives the following equations for 


{z) - exp {NSIp,p]) 
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W^{u,v) and U{^) 

w,iu,v) = pwpf; 


k=o ^ j dudvQ^i,{u,v\k){v/uY 


(29) 


[c/(m)]“^ = E 

where 


kpk f du dv Q^,{u, v\k — 1) {v/u )' 


^ j dudvQ^i,(u,v\k){v/u)'^^ 
Qti{u,v\k) = / n dUr dVr dJr W^(Ur, Vr) P{ Jr) 


X S 


\r—l 

1 


U — 




{z* + Er=l J?^r) 

m 


is the conditional distribution of u and v for a given de¬ 
gree k. Finally, we substitute eq. (l28l) in eq. (1271) and 
perform the limit n± —>■ ±^, from which the expression 
for the large N behavior of GNifJ-, A) is derived 

0Ar(/r,A)= lim expj - ^ln[C/(^)]^ I 
€^0+ I Z J 

{ oo r /* — 

ivE^feln / du dv Qf^{u, v\k) ^ 

k=o “ 

(31) 


In principle, eq. (IMI) determines completely the large- 
N behavior of the characteristic function for the index 
distribution of random graphs with arbitrary degree and 
edge distributions, as long as a solution for W^(u,r;) is 
extracted from the intricate self-consistent equation (l29l) . 

For A = 0, one can show that W^(u, v) = d{u—v)Rf_i{u) 
solves eq. (1^ . provided the normalized distribution 
R^(u) fulfills a certain equation, whose particular form 
is not relevant in this case. Thus, the characteristic func¬ 
tion at A = 0 simply reads 

0jv(m,O) = exp (32) 


which yields the delta peak Vn{K, 0) = 5 [K — iV/2] for 
the index distribution, after substituting eq. (15^ in eq. 
(0- This result reveals that, in order to access the index 
fluctuations in this case, one needs to compute the next- 
order contribution to Gn{p,0) for large N. The same 
situation arises in the replica approach for the GOE en¬ 
semble [ 3 ] . We present in the next section explicit results 
for the mean and the variance of the index for specific 
random graph models in the regime |A| > 0. 


VI. STATISTICAL PROPERTIES OF THE 
INDEX 

It is straightforward to check from eqs. ([9]) and (|^ 
that the moments (AT") scale as (AT") oc iV" for large N. 


In particular, the mean and the variance read 


(AT) =iVm(A), 

(K^) - (Kf = Na^{X), (33) 


where the prefactors m(A) and o-^(A) depend on the spe¬ 
cific graph ensemble via the distributions pk and P{J)- 
Equation (13311 differs strikingly from rotationally invari¬ 
ant ensembles of random matrices , where the 

variance of the typical index fluctuations is of O(lniV) 
and the prefactor is independent of A BIIMI. Erom 
eq. (15^ we conclude that cr^(O) = 0, which suggests 
that the index variance of random graphs with an arbi¬ 
trary degree distribution exhibits the logarithmic scaling 
(AT^) — (A')^ oc In V for large N at this particular A. This 
is confirmed below for the case of ER random graphs by 
means of numerical diagonalization results. 

For |A| > 0, the intensive quantities m(A) and cr^(A) 
are obtained directly from eqs. m and m, i.e., from the 
coefficients of the expansion of t/Ar(/i, A) around p = 0. 
In general, m(A) and cr^(A) are given in terms of averages 
with the distribution Wo(u, v) = lim^_>o IF^(m, v), whose 
self-consistent equation is derived by performing the limit 
^ 0 in eq. (l29)) 


tXJ J 

Wo{u,v) ='^—Qo{u,v\k - 1). (34) 

fc=o ^ 


The object Wo{u,v) may be interpreted as the averaged 
joint distribution of the diagonal resolvent elements at 
the two different points z and —z* of the complex plane. 
The resolvent elements at z and —z* are both calculated 
on the same cavity graph [IJ, , defined as the graph 
in which an arbitrary node and all its edges are deleted. 

Equation (I34p has a simpler form when compared to 
eq. and numerical solutions for Wo{u,v) can be 

obtained using the population dynamics algorithm [^ . 
where the distribution Wq{u,v) is parametrized by a 
large set {m^, containing M pairs of stochastic 

random variables. These are updated iteratively accord¬ 
ing to their joint distribution Wq{u,v), governed by eq. 
(IMl) . until Wo(m, v) attains a stationary profile. The limit 
e —>■ 0 + in eq. m is handled numerically by calculat¬ 
ing Wo{u, v) for small but finite values of e. We refer the 
reader to references [H, [IJ, [13 fo'' further details regard¬ 
ing the population dynamics algorithm in the context of 
random matrices and some technical points involved in 
the limit e —>■ O'*'. Since the eigenvalue distribution pat (A) 
is symmetric around A = 0, m(X) and cr^(A) obey the re¬ 
lations m(—A) = 1 — m{X) and cr^(A) = (t^(—A). Hence 
the results for m(A) and cr^(A) discussed below are lim¬ 
ited to the sector A > 0. 
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A. Erdos-Renyi random graphs 


For ER random graphs the quantities m(A) and cr^(A) 
read 


m(X) = lim 
£^• 0 + 


du dv du' dv' Wo(u, v) 


X Wo{u ,v')Ai{u,v,u\v') 


(35) 


cr^(A) = lim 
£->■ 0 + 


du dv du' dv' Wq{u, v) 


X Wo{u', v')A 2 {u,v,u', v') 


(36) 


where 


Ai(u, V] u', v') = ^ {Fj(u, V] u', v'))j “ ^ ’ 

A 2 (u,v\u', v') = {[Fj{u,v]u',v')f 


+ 


dTT^ 


/ v \ . 

/ v'\ 

1 

\ / vy 

-n- 


In (-) 

V u/ 

V / 


V uJ \ 


^ 2 


with 

Fj{u,v,u',v')=hY(j-FL^^. (37) 

The distribution Wq{u^v) is calculated numerically from 
eq. (1341) using the population dynamics algorithm with 

the degree distribution pk = " u!' of ER random graphs 

0 - 

In figures [T] and we present numerical results for 
to(A) and cr^(A) in the case of ER random graphs with 
P{J) = 6{J — 1). The discontinuous behavior of m{X) 
for small average degree c reflects the presence of delta 
peaks in the eigenvalue distribution, due to the proximity 
of the percolation transition (ssj. In fact, all connected 
components of ER random graphs are finite trees and the 
spectrum is purely discrete for c < 1, while the heights 
of these peaks decrease exponentially with increasing c 
[ 3 ^. The calculation of the integrated density of states 
presented here allows to determine, for A —>■ 00 , not 
only the location of the most important delta peaks in 
the spectrum, but also their relative weights, given by 
the size of the discontinuities of m(A). The exactness of 
our results for m(A) is confirmed by the comparison with 
numerical diagonalization data, as shown in figure [TJ 

The results for the prefactor (t^(A) of ER random 
graphs are shown in figure [21 Eor the smaller values 
of c, the index fluctuations are generally stronger and 
cr^(A) exhibits an irregular behavior, both features re¬ 
lated to strong sample to sample fluctuations of the 
graph structure close to the percolation critical point. 
The prominent feature of figure |2] is that o'^(A) shows 



Figure 1: Numerical results for the averaged intensive index 
m(A) of Erdos-Renyi random graphs with the distribution of 
edges P{J) = (5(J— 1), obtained using the population dynam¬ 
ics algorithm (solid lines) with M = 10® random variables and 
e = 10“®. Numerical diagonalization results (symbols), cal¬ 
culated from an ensemble of 100 matrices of size N = 3200, 
are shown as a comparison. 



Figure 2: Numerical results for the prefactor cr^(A) of the 
index variance of Erdos-Renyi random graphs with the distri¬ 
bution of edges P{J) = 5{J — 1), obtained using the popula¬ 
tion dynamics algorithm with M = 10® random variables and 
e = 10“®. 


a non-monotonic behavior, with a maximum for a cer¬ 
tain intermediate value of A and a vanishing behavior at 
A = 0, which signals the breakdown of the linear scaling 
(RT^) — (RT)^ oc N. This is confirmed by the numerical 
diagonalization results of figure [31 where (RT^) — {K)'^ is 
calculated as a function of N for c = 3. The results of 
figure [3(a) for different values of A > 0, display a linear 
behavior for increasing N, with slopes in full accordance 
with the theoretical values for (T^(A), as indicated on the 
caption. On the other hand, figure 3(b) shows that the 
index variance scales as (RT^) — (RT)^ oc hiN for A = 0, 
similarly to the behavior of rotationally invariant ensem¬ 
bles 
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(a)Index variance for A > 0. The solid lines represent the 
linear fit {K'^) — {K)^ = a + bN, with the values of the slope b 
indicated next to each straight line. The theoretical values for 
(t^(A), calculated through the numerical solution of eq. II36II . 
are given by cr^ (0.5) = 0.015, cr^ (3.0) = 0.0085 and 
fT2(3.5) = 0.0040. 



(b)Index variance for A = 0. The solid line represents the 
logarithmic fit {K'^) — (K)^ = a + blnN, with the slope 
b = 0.47(6). 

Figure 3: Numerical diagonalization results for the index vari¬ 
ance of Erdos-Renyi random graphs with c = 3 as a function 
of the number of nodes N. Each data point is calculated from 
an ensemble with S independent realizations of the adjacency 
matrix A, where S has been chosen according to S = . 

The solid lines represent the best fits of the numerical data. 


B. Random regular graphs 

In the case of random regula r g raphs, the degree dis¬ 
tribution is simply pk = 5k,c [17j |. where c > 2 is an 
integer. Firstly, let us consider the situation in which 
the values of the edges are fixed, i.e., their distribution 
reads P{J') = S{J' — J), with J G R. In this case, eq. 
(|2^ has the following solution for arbitrary p 




) 5{u-g)S{v + g*), (38) 

g*u; 


where g is a root of the algebraic equation 

(c-I)JV - + 1 = 0. (39) 


The quantity g represents the dia gon al elements of the 
resolvent on the cavity graph pll. l33|. Substituting eq. 
(1551) in eq. (1511) and using the above quadratic equation, 
we get 


0]\j(g,X)= lim exp[igNm(z)], 
£->0 + 


(40) 


where 

m{z) = —Im [in (z — cJ^g)] — [in (l — J^g^)] ■ 

(41) 

Equation (HUl) is the large-behavior of GN{g, A) for ran¬ 
dom regular graphs in the absence of edge fluctuations. 
By choosing the proper roots of eq. (155)) in the different 
sectors of the spectrum [s^, we can perform the limit 
lim(,_>g+ m(z) and derive the following analytical result 
for A > 0 


m(A) = 1-1— tan 

TT 


-c\/A^ - A^ 

A(c-2) 


-tan 

27r 


-1 


A\/Afa - A^ 

A2 -2c(c- 1)J2 


(42) 


with |Ab| = 2\J\y/c— 1 denoting the band edge of the 
continuous spectrum of random regular graphs [13, ■ 

Equation (1551) coincides with the average integrated den¬ 
sity of states in the bulk of a Cayley tree [s^ and it con¬ 
verges to the result for the GOE ensemble when c ^ 1 
[ 3 ], as long as we rescale J according to J —>■ J j y/c. 
The substitution of eq. OMl) in eq. © yields a delta 
peak Vn{K,X) = 5[K — Nm{X)], which implies that 
cr^(A) = 0. This suggests that the index variance ex¬ 
hibits the logarithmic scaling (/C^) — (A')^ oc IniV for 
arbitrary A. The latter property is consistent with the 
absence of localized states and the corresponding repul¬ 
sion between nearest-eigenvalues, which is common to the 
whole spectrum of random regular graphs with uniform 
edges [26142^ . 

The above results are clearly due to our trivial choice 
for P{J). The spectrum of random regular graphs con¬ 
tains localized states in the presence of edge disorder 
[13,113 and one can expect that cr^(A) exhibits a nontriv¬ 
ial behavior as long as P{J) has a finite variance. The 
functions ra[X) and cr^ (^) for random regular graphs with 
an arbitrary distribution P{J) read 


m(A) = lim \^Ki{z) - Li(z) 

ZTT £->0+ LZ 




+ [L,{z)f-L^{z)\, (43) 
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0 2 4 6 

A 

(b) 

Figure 4: Numerical results for the averaged intensive index 
m(A) and the prefactor cr^(A) of the index variance of random 
regular graphs with edges drawn from the Gaussian distribu¬ 
tion P{J) = (27r)~ 2 exp (— J^/2), obtained using the pop¬ 
ulation dynamics algorithm (solid lines) with M = 5 x 10® 
random variables and e = 10“®. Numerical diagonalization 
results (symbols), calculated from an ensemble of 100 matri¬ 
ces of size N = 4000, are shown as a comparison. 



(a)Erd6s-Renyi random graphs with the distribution of the 
edges P(J) = S(J — 1). 



(b)Regular random graphs with the distribution of the edges 
P(J) = (27r)“2 exp [—J^l2). 

Figure 5: Numerical diagonalization results (symbols) for the 
distribution of the intensive index of random graphs with 
c = 5 and A = 1. The histograms were generated from 10® 
independent samples for the intensive index of the adjacency 
matrix A. The solid lines are Gaussian distributions with 
mean and variance taken from the data. 


where Kn{z) and L„(z) are calculated from 

Kn{z) = J dudvdu'dv'Qo{u,v\c—l)Qo{u\v'\c— 1) 
X {[Fj{u,v,u',v’)y")j , 


Ln{z) = J dudvQo{u,v\c) In ^-^ 


(44) 


Figure |4] shows population dynamics results for m(A) and 
(t^(A) in the case of a Gaussian distribution P{J) = 
(27r)~^ exp (—J^/2). The function m(A) does not dis¬ 
play any noticeable discontinuity, as observed previously 
for ER random graphs, due to the absence of discon¬ 
nected clusters in the case of large random regular graphs 
17|. In addition, we note that cr^(A) has qualitatively the 
same non-monotonic behavior as in ER random graphs, 
exhibiting a maximum for a certain A and approaching 
zero as A —> 0. Numerical diagonalization results for 
large matrices A, also shown in figured confirm the cor¬ 


rectness of our theoretical approach. 


C. The index distribution 

In this subsection, we inspect the full index distribu¬ 
tion of random graphs using numerical diagonalization, 
instead of undertaking the more difficult task of calculat¬ 
ing the characteristic function from the numerical solu¬ 
tion of eqs. (1^ and (IHTl) . We restrict ourselves to A > 0, 
where the index variance scales linearly with N ^ 1. 

In figure[5]we show results for the distribution pAr(/c, A) 
of the intensive index = /Civ(A)/7V in the case of 

ER and random regular graphs with c = 5, obtained from 
numerical diagonalization for A = 1. For each value of 
N, the results are compared with a Gaussian distribution 
(solid lines) with mean and variance taken from the data, 
which confirms the Gaussian character of the typical in¬ 
dex fluctuations for both random graph models when N 
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is large but finite. 

Overall, our results suggest that, for 1 and |A| > 
0, the intensive index of ER and random regular graphs 
is distributed according to 


N \ N 2 

( 45 ) 

with non-universal parameters cr^(A) and m(A) that de¬ 
pend on the underlying random graph model as well as 
on the particular value of the threshold A. The func¬ 
tion pN{k,X) converges to pN{k,X) = i5[A: — m(A)] for 
N ^ oo, but the rate of convergence is slower when com¬ 
pared to rotationally invariant ensembles , due 

to the logarithmic scaling of the index variance with re¬ 
spect to N in the latter case. On the other hand, the 
Gaussian nature of the index fluctuations for ^ 1 
seems to be an ubiquitous feature of random matrix mod¬ 
els. 


VII. FINAL REMARKS 

We have presented an analytical expression for the 
characteristic function of the index distribution describ¬ 
ing a broad class of random graph models, which com¬ 
prises graphs with arbitrary degree and edge distribu¬ 
tions. Ideally, this general result gives access to all mo¬ 
ments of the index distribution in the limit V —>■ oo. We 
have shown that the index variance of typical fluctuations 
is generally of 0{N), with a prefactor fT^(A) that depends 
on the random graph model under study as well as on the 
threshold A that defines the index through eq. Q. In 
particular, cr^(A) follows an intriguing non-monotonic be¬ 
havior for random graphs with localized eigenstates: it 
exhibits a maximum at a certain |A| > 0 and a vanish¬ 
ing behavior at A = 0. Numerical diagonalization data 
confirm the theoretical results and support the Gaussian 
form of the typical index distribution for the random 
graphs considered here (see eq. (H5|) b completing the 
picture about the index statistics. 

Our results differ with those of rotationally invariant 
ensembles, where the index variance is of 0(ln7V), with 
a prefactor that is independent of A and has an univer¬ 
sal character. We argue that this difference in the scal¬ 
ing forms arises due to the presence of localized states 
in the spectrum of some random graphs. In the local¬ 
ized sectors, the eigenvalues do not repel each other and 
behave as uncorrelated random variables, such that the 
total number of eigenvalues contained in finite regions 
within the localized phase suffers from stronger finite size 
fluctuations as compared to regions within the extended 


phase, where level-repulsion tends to equalize the space 
between neighboring eigenvalues. On the other hand, the 
Gaussian nature of typical index fluctuations seems to be 
a robust feature of random matrix models. 

On the methodological side, the replica approach as 
devised here departs from the representation of the char¬ 
acteristic function in terms of real Gaussian integrals, 
instead of the fermionic Gaussian integrals adopted in 
reference Q- In the situations where cr^(A) = 0, the log¬ 
arithmic scaling of the index variance is obtained in our 
setting from the next-to-leading order terms, for large 
V, in the saddle-point integral of eq. (1^ . These contri¬ 
butions come from 0(1/VN) fluctuations of the order- 
parameter and they are handled following the ideas of 
reference [s^. Indeed, we have precisely recovered the 
analytical results for the GOE ensemble Q employing 
this strategy [i^ , and the same approach can be used to 
calculate the prefactors in situations where the variance 
of random graphs is of 0(ln V). 

Our work opens several perspectives in the study of the 
typical index fluctuations. Firstly, it would be worth hav¬ 
ing approximate schemes or numerical methods to solve 
eq. (EH) and obtain the distribution Wfj,{u,v), which 
would allow to fully determine the characteristic func¬ 
tion for random graphs. Due to the versatile character 
of the replica method, the study of the averaged inte¬ 
grated density of states of the Anderson model on regu¬ 
lar graphs and its sample to sample fluctuations is 
just around the corner. It would be also interesting to 
inspect the robustness of the Gaussian form of the index 
fluctuations in random matrix ensembles with strong in¬ 
herent fluctuations, such as Levy random matrices 
and scale-free random networks [1^. The index statis¬ 
tics of both random matrix models can be treated using 
the replica approach as developed here. In fact, scale-free 
random graphs, crucial in modelling many real-world net¬ 
works appearing in nature [1^, can be studied directly 
from our work by choosing the degree distribution as 
Pk ~ k~'^ (2 < 7 < 3), which yields random graphs with 
strong sample to sample degree fluctuations. Finally, we 
point out that the different scaling behaviors of the in¬ 
dex variance should have important consequences to the 
relaxation properties and search algorithms on complex 
energy surfaces. 
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